library(dplyr) # makes variables
library(readr) # reads in dataset
library(broom) # extracts statistics from models
library(ggplot2) # makes plots
library(ggpubr) # arranges plots
library(BayesFactor) # calculates Bayes Factors
library(psych) # calculates reliability statistics
library(stargazer) # makes tables

# Set this to the location of the dataset on your machine
setwd("~/Data/COVID-19_Accountability")

#### IMPORT DATA ####

d <- read_csv("covid_serious_Moniz_2020.csv")

#### VARIABLE CONSTRUCTION ####

# condition indicator constructed as follows:
d <- mutate(d, condition = case_when(
      Manipulations_DO == "control" ~ 0,
      Manipulations_DO == "minimize" ~ 1,
      Manipulations_DO == "minimize.context" ~ 2,
      Manipulations_DO == "amplify" ~ 3
))

# R says he/she had heard of treatment before
d$pretreated.binary <- ifelse(d$pretreated <= 2, 1, 0)

# Problem seriousness DV (scaled 0-1), larger values = deaths are a bigger problem
d$deaths.bigproblem <- ((6 - d$deaths.problem) - 1) / 4

# Trump Job approval scale DV
d <- mutate(d, trump.job = (as.numeric(d$trump.weak) + as.numeric(d$trump.incomp) +
                                  as.numeric(d$trump.effective)
                            + as.numeric(d$trump.helpful)) / 4)
d <- mutate(d, trump.job = 10 - trump.job) # higher = better job
d <- mutate(d, trump.job = (trump.job - 1)/ 8) # scale from 0-1

# policy support scale DV; scaled 0-1; larger = more support 
d <- mutate(d, policy.scale =( ((8 -((Q47 + Q60 + Q61 + Q62 + Q63 + Q64 + Q65 + Q66)/8))) - 1) /6)

# risky health-behavior intentions DV
d$face.mask <- (8 - d$Q75) # reverse-coded
# larger values = more risky intentions
d <- mutate(d, intentions =( (((Q67+ Q77+ Q76+ Q72+ face.mask+ Q73+ Q74+ Q69)/8)) - 1) /6)


#### TEST RELIABILITY OF SCALES ####


trump.job <- data.frame(weak = as.numeric(d$trump.weak), incomp = as.numeric(d$trump.incomp),
                        effec = as.numeric(d$trump.effective),
                        help = as.numeric(d$trump.helpful))
psych::alpha(trump.job)

policy <- d %>% select(Q47, Q60, Q61, Q62, Q63, Q64, Q65, Q66)
psych::alpha(policy)

intentions <- d %>% dplyr::select(Q67, Q77, Q76, Q72, face.mask, Q73, Q74, Q69)
psych::alpha(intentions)



####################
#### INFERENCES ####
####################

#### HYPOTHESIS ONE ####
## DIFFERENCES IN PROBLEM STATUS ##
summary(lm(deaths.bigproblem ~ factor(condition), d))

t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 1, ], 
       alternative = "greater", conf.level = 1-(.05/4))

t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 3, ], 
       alternative = "less", conf.level = 1-(.05/4))


# Bayes Factors
# nullInterval argument specifies the one-sided hypothesis test 
# (difference can be between 0 and infinity, not negative)
ttestBF(d$deaths.bigproblem[d$condition == 0], d$deaths.bigproblem[d$condition == 1],
        nullInterval = c(0, Inf))
ttestBF(d$deaths.bigproblem[d$condition == 0], d$deaths.bigproblem[d$condition == 3],
        nullInterval = c(0, Inf))



#### Figure 1 ####
## Contains all four DVs ##

## Perceived Seriousness of Death Toll DV
# estimating standard errors of means for CIs
mean0 <- broom::tidy(t.test(d[d$condition == 0,"deaths.bigproblem"], conf.level = .9875))

mean1 <- broom::tidy(t.test(d[d$condition == 1,"deaths.bigproblem"], conf.level = .9875))

mean2 <- broom::tidy(t.test(d[d$condition == 2,"deaths.bigproblem"], conf.level = .9875))

mean3 <- broom::tidy(t.test(d[d$condition == 3,"deaths.bigproblem"], conf.level = .9875))

h1means <- bind_rows(mean0, mean1, mean2, mean3)

h1means$group <- factor(c("Control", "Minimize", "Min. + Counter-Info", "Amplify"),
                        levels = c("Control", "Minimize", 
                                   "Amplify", "Min. + Counter-Info"))


serious.plot <- h1means %>% ggplot(aes(group, estimate)) + geom_col(width = .35) + 
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
      xlab("") + ylab("Mean Response") + 
      theme_minimal() + ylim(c(0, .8)) + 
      labs(subtitle = "Panel A: Perceived Seriousness of COVID-19 Death Toll") +
      theme(plot.subtitle = element_text(hjust = 0.5))


## Trump job performance DV 
# estimating standard errors of means for CIs
mean0 <- broom::tidy(t.test(d[d$condition == 0,"trump.job"], conf.level = .9875))

mean1 <- broom::tidy(t.test(d[d$condition == 1,"trump.job"], conf.level = .9875))

mean2 <- broom::tidy(t.test(d[d$condition == 2,"trump.job"], conf.level = .9875))

mean3 <- broom::tidy(t.test(d[d$condition == 3,"trump.job"], conf.level = .9875))

h2means <- bind_rows(mean0, mean1, mean2, mean3)

h2means$group <- factor(c("Control", "Minimize", "Min. + Counter-Info", "Amplify"),
                        levels = c("Control", "Minimize", 
                                   "Amplify", "Min. + Counter-Info"))

trump.job.plot <- h2means %>% ggplot(aes(group, estimate)) + geom_col(width = .35) + 
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
      xlab("") + ylab("") + 
      theme_minimal() + ylim(c(0, .8)) +
      labs(subtitle = "Panel B: Trump Performance Evaluation") +
      theme(plot.subtitle = element_text(hjust = 0.5))

## Support for Social-Distancing Policies DV
# estimating standard errors of means for CIs
mean0 <- broom::tidy(t.test(d[d$condition == 0,"policy.scale"], conf.level = .9875))

mean1 <- broom::tidy(t.test(d[d$condition == 1,"policy.scale"], conf.level = .9875))

mean2 <- broom::tidy(t.test(d[d$condition == 2,"policy.scale"], conf.level = .9875))

mean3 <- broom::tidy(t.test(d[d$condition == 3,"policy.scale"], conf.level = .9875))

h3means <- bind_rows(mean0, mean1, mean2, mean3)

h3means$group <- factor(c("Control", "Minimize", "Min. + Counter-Info", "Amplify"),
                        levels = c("Control", "Minimize", 
                                   "Amplify", "Min. + Counter-Info"))

policy.plot <- h3means %>% ggplot(aes(group, estimate)) + geom_col(width = .35) + 
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
      xlab("") +  ylab("") + 
      theme_minimal() + ylim(c(0, .8)) +
      labs(subtitle = "Panel C: Social Distancing Policy Support") +
      theme(plot.subtitle = element_text(hjust = 0.5))

## Intent toward Risky Health Behaviors DV
# estimating standard errors of means for CIs
mean0 <- broom::tidy(t.test(d[d$condition == 0,"intentions"], conf.level = .9875))

mean1 <- broom::tidy(t.test(d[d$condition == 1,"intentions"], conf.level = .9875))

mean2 <- broom::tidy(t.test(d[d$condition == 2,"intentions"], conf.level = .9875))

mean3 <- broom::tidy(t.test(d[d$condition == 3,"intentions"], conf.level = .9875))

h4means <- bind_rows(mean0, mean1, mean2, mean3)

h4means$group <- factor(c("Control", "Minimize", "Min. + Counter-Info", "Amplify"),
                        levels = c("Control", "Minimize", 
                                   "Amplify", "Min. + Counter-Info"))

intentions.plot <- h4means %>% ggplot(aes(group, estimate)) + 
      geom_col(width = .35) + 
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
      xlab("") + ylab("") + 
      theme_minimal() + ylim(c(0, .8)) +
      labs(subtitle = "Panel D: Risky Health Behavior Intentions") +
      theme(plot.subtitle = element_text(hjust = 0.5))

fig1 <- ggarrange(serious.plot, trump.job.plot, policy.plot, intentions.plot,
                  widths = c(2))

annotate_figure(fig1,
                bottom = text_grob("Brackets represent 98.75% confidence intervals.",
                                         size = 8))
ggsave("fig1.png", scale = .66)


#### HYPOTHESIS TWO ####
## DIFFERENCES IN TRUMP APPROVAL ##

with(d[d$condition == 0 | d$condition == 1,], t.test(trump.job ~ condition, 
                                                     alternative = "less", 
                                                     conf.level = 1-(.05/4)))
# Group sizes to calculate Bayes Factor
sum(!is.na(d$trump.job[d$condition == 0]))
sum(!is.na(d$trump.job[d$condition == 1]))

#Bayes Factor
ttest.tstat(-1.8413, 421, 405, simple = T, nullInterval = c(-Inf, 0))




#### HYPOTHESIS THREE ####
## DIFFERENCES IN TRUMP APPROVAL ##
## Between Problem-minimizing and Problem-minimizing with Counterargument ##

with(d[d$condition == 0 | d$condition == 2,], t.test(trump.job ~ condition, 
                                                     alternative = "less", 
                                                     conf.level = 1-(.05/4)))
with(d[d$condition == 1 | d$condition == 2,], t.test(trump.job ~ condition, 
                                                     alternative = "less", 
                                                     conf.level = 1-(.05/4)))

summary(lm(trump.job ~ factor(condition), d))


#### HYPOTHESIS FOUR ####

## DIFFERENCES IN POLICY ATTITUDES ##

with(d[d$condition == 0 | d$condition == 1,], t.test(policy.scale ~ condition,
                                                     alternative = "greater", 
                                                     conf.level = 1-(.05/4)))
with(d[d$condition == 0 | d$condition == 3,], t.test(policy.scale ~ condition,
                                                     alternative = "less", 
                                                     conf.level = 1-(.05/4)))

#Bayes Factors
# Group sizes are calculated thus and placed directly into BF code
sum(!is.na(d$policy.scale[d$condition == 0]))

ttest.tstat(-.31193, sum(!is.na(d$policy.scale[d$condition == 0])),
            sum(!is.na(d$policy.scale[d$condition == 1])), 
            simple = T, nullInterval = c(-Inf, 0))
ttest.tstat(.6366, sum(!is.na(d$policy.scale[d$condition == 0])),
            sum(!is.na(d$policy.scale[d$condition == 3])), 
            simple = T, nullInterval = c(0, Inf))



## DIFFERENCES IN RISK INTENTIONS ##

with(d[d$condition == 0 | d$condition == 1,], t.test(intentions ~ condition,
                                                     alternative = "greater", 
                                                     conf.level = 1-(.05/4)))
with(d[d$condition == 0 | d$condition == 3,], t.test(intentions ~ condition,
                                                     alternative = "less", 
                                                     conf.level = 1-(.05/4)))

#Bayes Factors
# Group sizes are calculated directly in BF code
ttest.tstat(-.15943, sum(!is.na(d$intentions[d$condition == 0])),
            sum(!is.na(d$intentions[d$condition == 1])), 
            simple = T, nullInterval = c(-Inf, 0))
ttest.tstat(.005175, sum(!is.na(d$intentions[d$condition == 0])),
            sum(!is.na(d$intentions[d$condition == 3])), 
            simple = T, nullInterval = c(0, Inf))


#### EXPLORATORY ANALYSES ####

## DIFFERENCES IN PROBLEM STATUS ##

t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 2, ])

# Bayes Factors
# Group sizes are calculated directly in BF code
ttest.tstat(4.1146, sum(!is.na(d$deaths.bigproblem[d$condition == 0])),
            sum(!is.na(d$deaths.bigproblem[d$condition == 2])), 
            simple = T)

t.test(trump.job ~ condition, d[d$condition == 0 | d$condition == 2, ], 
       conf.level = 1-(.05/4))



### MEDIATION ANALYSES ###

## HYPOTHESIS TWO ##
# Model: Perceived seriousness mediates treatment effect on Trump's job eval.

# Using PROCESS
# Problem-Minimizing condition
process(data = d[d$condition == 0 | d$condition == 1,], 
        y = "trump.job", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"), model = 4, conf = 98.75, seed = 155292)

# Problem-Amplifying condition
process(data = d[d$condition == 0 | d$condition == 3,], 
        y = "trump.job", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"), model = 4, conf = 98.75, seed = 155292)


## HYPOTHESIS THREE ##

process(data = d[d$condition == 0 | d$condition == 2,], 
        y = "trump.job", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"), model = 4, conf = 98.75, seed = 155292)

# Sensitivity Analysis Using Imai et al.'s Causal Mediation Analysis package
# This package requires equal-size datasets be used for outcome and mediator models, so
# 1 observation with missing value on the DV needs to be dropped manually.
library(mediation)
med.data <- d[!is.na(d$trump.job),]
# Linear model of mediator
modm <- lm(deaths.bigproblem ~ condition + ideology + follow.news, 
           data = med.data[med.data$condition == 0 | med.data$condition == 2,])
# Linear model of outcome
mody <- lm(trump.job ~ deaths.bigproblem + condition + ideology + follow.news, 
           data = d[d$condition == 0 | d$condition == 2,])
# Estimates the mediation model
cma <- mediate(modm, mody, treat = "condition", mediator = "deaths.bigproblem",
                covariates = c("ideology", "follow.news"))
# Conducts sensitivity analysis
sencma <- medsens(cma)
summary(sencma)

# Calculating sample correlations for comparing sensitivity to observed correlations
cor.test(d$deaths.bigproblem, d$ideology)
cor.test(d$deaths.bigproblem, d$follow.news)
cor.test(d$trump.job, d$ideology)
cor.test(d$trump.job, d$follow.news)

## HYPOTHESIS FOUR ##

# DV = policy support
process(data = d[d$condition == 0 | d$condition == 1,], 
        y = "policy.scale", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"),
        model = 4, conf = 98.75, seed = 155292)

process(data = d[d$condition == 0 | d$condition == 3,], 
        y = "policy.scale", x = "condition", m = "deaths.bigproblem", 
        cov = c("ideology", "follow.news"),
        model = 4, conf = 98.75, seed = 155292)

# DV = behavioral intentions
process(data = d[d$condition == 0 | d$condition == 1,], 
        y = "intentions", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"),
        model = 4, conf = 98.75, seed = 155292)

process(data = d[d$condition == 0 | d$condition == 3,], 
        y = "intentions", x = "condition", m = "deaths.bigproblem",
        cov = c("ideology", "follow.news"),
        model = 4, conf = 98.75, seed = 155292)


## PRETREATMENT ESTIMATE ##

# 1 = seen/heard of stimulus once or several times
prop.table(table(d$pretreated.binary[d$condition == 0]))



#### APPENDIX TABLES ####

# TABLE 1: Descriptive statistics by condition
# Calculate variable means by condition
means <- d %>% group_by(condition) %>%
      summarize(ageM = mean(age, na.rm=T),
                        ideologyM = mean(ideology, na.rm=T),
                        follow.newsM = mean(follow.news, na.rm = T))

# Calculate variable standard deviations by condition
sds <- d %>% group_by(condition) %>%
                summarize(ageSD = sd(age, na.rm=T),
                ideologySD = sd(ideology, na.rm=T),
                follow.newsSD = sd(follow.news, na.rm=T))

# round to 2 decimal places
sds <- round(sds, 2)
means <- round(means, 2)

# write in names of conditions
means[,1] <- c("Control", "Minimize", "Minimize + Counter-Info", "Amplify")
# give columns names
colnames(means) <- c("Condition", "Age", "Ideology", "Pol. Interest")

# print table of means (added SDs by hand)
stargazer(means, summary = F) # table output


# TABLE 2: Calculating means, difference of means confidence intervals, 
# and Ns for experimental groups
# These were then entered by hand into a LaTeX table

# Calculate means and difference-in-means standard errors for problem seriousness DV
t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 1, ], 
alternative = "greater", conf.level = 1-(.05/4))
t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 2, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(deaths.bigproblem ~ condition, d[d$condition == 0 | d$condition == 3, ], 
       alternative = "greater", conf.level = 1-(.05/4))
# Calculate Ns for problem seriousness DV
sum(!is.na(d$deaths.bigproblem[d$condition == 0]))
sum(!is.na(d$deaths.bigproblem[d$condition == 1]))
sum(!is.na(d$deaths.bigproblem[d$condition == 2]))
sum(!is.na(d$deaths.bigproblem[d$condition == 3]))

# Calculate means and difference-in-means standard errors for Trump performance DV
t.test(trump.job ~ condition, d[d$condition == 0 | d$condition == 1, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(trump.job ~ condition, d[d$condition == 0 | d$condition == 2, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(trump.job ~ condition, d[d$condition == 0 | d$condition == 3, ], 
       alternative = "greater", conf.level = 1-(.05/4))
# Calculate Ns for Trump performance DV
sum(!is.na(d$trump.job[d$condition == 0]))
sum(!is.na(d$trump.job[d$condition == 1]))
sum(!is.na(d$trump.job[d$condition == 2]))
sum(!is.na(d$trump.job[d$condition == 3]))

# Calculate means and difference-in-means standard errors for policy support DV
t.test(policy.scale ~ condition, d[d$condition == 0 | d$condition == 1, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(policy.scale ~ condition, d[d$condition == 0 | d$condition == 2, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(policy.scale ~ condition, d[d$condition == 0 | d$condition == 3, ], 
       alternative = "greater", conf.level = 1-(.05/4))
# Calculate Ns for Trump policy support DV
sum(!is.na(d$policy.scale[d$condition == 0]))
sum(!is.na(d$policy.scale[d$condition == 1]))
sum(!is.na(d$policy.scale[d$condition == 2]))
sum(!is.na(d$policy.scale[d$condition == 3]))

# Calculate means and difference-in-means standard errors for Trump behavioral intentions DV
t.test(intentions ~ condition, d[d$condition == 0 | d$condition == 1, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(intentions ~ condition, d[d$condition == 0 | d$condition == 2, ], 
       alternative = "greater", conf.level = 1-(.05/4))
t.test(intentions ~ condition, d[d$condition == 0 | d$condition == 3, ], 
       alternative = "greater", conf.level = 1-(.05/4))
# Calculate Ns for Trump behavioral intentions DV
sum(!is.na(d$intentions[d$condition == 0]))
sum(!is.na(d$intentions[d$condition == 1]))
sum(!is.na(d$intentions[d$condition == 2]))
sum(!is.na(d$intentions[d$condition == 3]))


## TABLE 3: Trump pandemic-specific performance evaluation
# recode variables from 0 to 1
trump.job <- data.frame(weak = as.numeric((9 - d$trump.weak)/8), 
               incomp = as.numeric((9 - d$trump.incomp)/8),
                        effec = as.numeric((9 - d$trump.effective)/8),
                        help = as.numeric((9 - d$trump.helpful)/8))
# print table
stargazer(trump.job, summary = T, digits = 2, summary.stat = c("mean", "sd"))


## TABLE 4: Social-distancing policy support
# # recode variables from 0 to 1
policy <- mutate(policy, (8 - !!all_vars(policy))/7 )
policy <- data.frame(policy)
# print table
stargazer(policy, summary = T, digits = 2, summary.stat = c("mean", "sd"))


## TABLE 5: risky health-behavior intentions
# # recode variables from 0 to 1
intentions <- mutate(intentions, (!!all_vars(intentions))/7)
intentions <- data.frame(intentions)
# print table
stargazer(intentions, summary = T, digits = 2, summary.stat = c("mean", "sd"))

## TABLE 6: experimental stimuli
# Text is taken from the questionnaire

